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I. INTRODUCTION 

The problem of heat transfer, electronic, phononic 
and photonic, in molecules and nanosystems has recently 
gained lots of interest [1-7]. In molecules, understanding 
heat flow is crucial for controlling reactivity, molecular 
dynamics, and kinetics [8]. In nanosystems, heat transfer 
has recently attracted much attention with implications 
in thermal machinery [9-11], information processing and 
computation [12, 13], and molecular-based thermoelec- 
tricity [14-16]. Of special interest are hybrid structures, 
e.g. normal metal-superconductor junctions with appli- 
cations in thermometry and refrigeration [17], and atom- 
radiation field systems, serving as a prototype for study- 
ing thermodynamics of quantum systems [18, 19]. 

From the theoretical point of view systems of interest 
include collections of bosons, fermions, spins, and mixed- 
statistics models [20]. For example, heat transfer from a 
dielectric solid into a molecule may be studied using a 
spin-boson model where the molecule is represented by a 
single anharmonic mode (spin) and the bulk includes a 
collection of harmonic modes (boson) [21]. In the analo- 
gous spin-fermion model an electronic excitation is trans- 
ferred between two metals through a local mode, model- 
ing a vibrating molecule. If the central mode is harmonic, 
the model may further describe radiative heat transfer 
between electronic conductors [7, 22] . 

In order to perform first principle quantum-mechanical 
calculations of heat transfer in nanosystems it is neces- 
sary to consider a model-independent non-perturbative 
definition of the heat current. This expression should be 
applicable in non-stationary cases, as well as in steady- 
state situations. While there is no unique definition of 
the heat current operator in non-relativistic systems [23] , 
the constructed expression should still fulfill a symme- 
try requirement, as we discuss below. We present here 
a consistent definition for the heat flux operator using 
a generic one-dimensional (ID) Hamiltonian. We show 
that this expression is useful for studying vibrational, 
electronic and spin mediated heat transfer, and that it 
yields a non-perturbative expression for the heat current 
in hybrid systems, e.g. at a solid-molecule-solid interface 



represented by a two-bath spin-boson model. 

Furthermore, the definition also brings in a useful phys- 
ical insight: We derive a necessary condition for energy 
conservation in various systems, bosonic and electronic, 
by calculating the commutator of the total flux opera- 
tor with H , the total Hamiltonian. If the current is a 
conserved quantity, the transport is ballistic, the con- 
ductivity diverges, and Fourier's law of heat conduction 
cannot he fulfilled [24]. 

Derivation of the Fourier's law from fundamental prin- 
ciples, classical [24-27], or quantum [28-30], is a great 
challenge in theoretical physics. Model calculations man- 
ifested that the onset of diffusional behavior delicately 
depends on the details of the system. It is still not clear 
what necessary and sufficient conditions must the Hamil- 
tonian fulfill for showing the Fourier's dynamics. Here we 
circumvent this challenge, and rather than test the ap- 
plicability of the Fourier's law in specific systems, derive 
a general, necessary condition for current conservation. 
Systems that do not obey this condition may satisfy the 
Fourier's law. As an example, we verify that in systems of 
harmonic oscillators the total heat current is conserved, 
so that once prepared, a current in a closed loop system 
will never vanish. 

Another implication of the proper definition of the heat 
current is the identification of a microscopic expression 
for the thermal conductivity in terms of Hamiltonian pa- 
rameters. This expression might be useful for studying 
the thermal conduction properties of molecular wires and 
spin chains. 

The paper is organized as follows. In Section II we 
discuss the general definition of the heat flux operator 
in one dimension. Section III applies this expression to 
complex structures, e.g. the spin-boson model and the 
spin-fermion model, prototype models for studying heat 
transfer in hybrid systems. In Section IV we show that 
a current conservation condition naturally emerges from 
the heat flux definition for both bosonic and fermionic 
Hamiltonians. Section V further explores current conser- 
vation in general ID systems. From the heat flux expres- 
sion the discrete Fourier's law can be naturally identified, 
as shown in Section VI. In Section VII we conclude. 
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II. DEFINITION OF THE ENERGY FLUX 



where 



OPERATOR FOR A GENERAL HAMILTONIAN 
IN ONE DIMENSION 



jfi,+i = |p°-/i°+i,ns,s+i)], (6) 



Defining a heat flux operator for a specific system 
such as phonons dates back to Hardy's early work [31]. 
The idea was appUcd to spin chains, see e.g. [32-35], 
and other ID systems, see e.g. Ref. [29]. A general 
flux (current) operator may be obtained by assuming 
that there exists an operator continuity equation, for 



instance 



^ dh(x,t) ^ dj{x,t) 



at 



dx 



in one dimension, where 
h{x,t) is the energy density operator and j{x,t) is the 
heat flux operator. For an A'^-site chain with M-states 
at each site, one can introduce a workable definition of 
the energy density operator, h{x,t) = X]g/is<5(a; — Xg), 
where hs is a discrete energy density operator of the 
*th site. The total Hamiltonian of the chain is therefore 
given hy H = J dxh{x,t) = X^g/is- Similarly, the heat 
flux can be written as the a of localized contributions 
j{x, t) = is5{x — Xg), so that the continuity equation 
can be written in a discrete form, 



dhg 
~dt 



(1) 



where a is the lattice spacing and jg/a is the current op- 
erator. The time evolution of hg in the Heisenberg rep- 
resentation satisfies the Heisenberg equation of motion, 
hg = i[H, hg], assuming hg does not depend on time in 



dt 



the Schrodinger representation {%= 1). This yields 



dhg -v^r, 

k 



(2) 



In general, Eq. (2) cannot be expressed in terms of the 
difference of two operators at two sites as in (1), yet we 
can identify the currents jg and jg-i for a specific model 
Hamiltonian. We use here a generic ID Hamiltonian with 
up to two-body nearest-neighbor interactions, 



N 



(3) 



where /i^ is the local Hamiltonian at site s. While the 
local energy density can not be imiquely defined [23], 
one could make a reasonable separation of V, and assign 
mixed terms half to site s and half to s -|- 1. With this 
partition the energy density at the sth local site becomes 



hg 



1 



K + T^[V{s,s + l) + V{s-l,s)]. 



(4) 



ihg, as required, when 



This equation satisfies H = ^ 

one sets V{N,N + 1) = 1/(0, 1) = 0. For this Hamilto- 
nian, the heat flux operator can be identified as 



7 - l^^) + l(^) 
Js — Js— >s+l ' Js ' 



(5) 



is a two-site contribution and 



ir = -^{[V{s,s+l\V{s+l,s + 2)] 
+ [V{s-l,s),V{s,s + l)]} 



(7) 



is an operator connecting four sites, accounting for higher 

order inter-site interaction terms. As we show below, 
in some cases it is exactly zero. It is also noticeable 
that in our case Eq. (2) could be written in terms of 
the difference of operators at two neighbor sites. The 
definition also naturally classifies the perturbative orders 
with respect to the inter-site coupling V: The order of 
the fiux operator (7) is higher than that of (6). 

The definitions (5)- (7) possess significant symmetric 
features. For instance, jg_^g_^.i trivially shows the ex- 
change symmetry jg^^^^ = — j^^^^g, assuming V{s, s + 
1) = V{s + 1, s). The exchange symmetry is an essential 
requirement when defining a current operator, since the 
current in opposite directions must have the same ab- 
solute value. The operators ji*^ has a similar exchange 
property, but four sites are involved. 

The definition (5) is state- and symmetry-independent 
unlike the expression utilized in Refs. [29, 36-38], jg = 
ia[h%V{s,s + 1)], which requires that the Hamiltonian 
fulfills the symmetric condition [/),2 + /j-^+i, V{s, s + l)] = 
0, see Appendix A for details. In order to increase gen- 
erality, Ref. [36] further suggests a 'symmetrized local 
flux' that has the same form as j^^Xg+\ ■ 

The heat flux operator was also extensively examined 
in ID chains in the absence of an on-site energy term 
{hPg = 0), e.g. the Heisenberg model at zero magnetic 
field [39] . In this case the energy at each site was defined 
as hg = V{s, s+l), leading to the current operator jg-i = 
ia[V{s — l,s),V{s,s + 1)]. Since in this paper we are 
interested in the opposite limit, i.e. in structures where 
the inter-site interaction is considered as a perturbation 
to the local energy, e.g. impurity models, the choice (4) 
for the local energy is more appropriate. 

Note, that we could have also defined a high order local 
interaction term U{s). For phononic systems U includes 
on-site interactions, incorporating harmonic and anhar- 
monic contributions. For fermionic systems U may rep- 
resent a local electron-electron repulsion. The potentials 
V{s, s+l) and V{s, s+l)+U{s)+U{s+l) indeed produce 
different flux operators. We adopt here the convention 
that local s interactions (one-body and many-body) are 
all included within the potential h'^. 

Finally, one could consider next nearest neighbor in- 
teractions, and by following the same procedure, identify 
the current jg . 
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III. CURRENT OPERATOR IN HYBRID 
STRUCTURES 

The definition (5) can be applied to non-identical in- 
teracting systems which are spaciously connected. For 
example, we may consider an impurity spin coupled to 
two solids, and study the heat current at the contact. 
The bulk, serving as a thermal reservoir, may be com- 
posed of electrons (the Kondo problem) [40], collections 
of harmonic modes (the spin-boson model) [41]. or spins 
[42, 43]. This impurity-bath scenario is the standard in 
molecular electronics and nanomcchanical experiments, 
where the heat transfer properties of a molecule con- 
nected to solid or liquid interfaces are investigated [2-5]. 

The generic impurity-bath Hamiltonian includes a cen- 
tral unit Hspin, two independent reservoirs H° {v = L, R) 
maintained at different temperatures, and system-bath 
couplings Vj/. The heat flux operator, e.g. at the L con- 
tact is given by Eqs. (5)-(7), disregarding for convenience 
the lattice constant a. Assuming that [Vl,Vr] = we 
find that the current from the L contact into the junc- 
tion is given by 

jL = ^[Hl-H,pin,VL]. (8) 

We apply next this result assuming either bosonic baths 
or electronic reservoirs. 

Spin-Boson model. — A two-level system connected to 
two harmonic baths held at different temperatures serves 
as a prototype model for investigating phononic transfer 
in a nonlinear molecular junction. Calculations at the 
level of the Master equation, assuming weak system-bath 
couplings while ignoring coherence effects, have revealed 
interesting dynamics, e.g. thermal rectification [21], neg- 
ative differential resistance [44], and pumping of heat 
[45]. It is of interest to derive a general expression for 
the heat current which is not limited to the weak cou- 
pling limit. Such an expression will open the door for 
non-perturbative calculations of heat current in strongly 
coupled molecular systems. The multi-bath spin-boson 
Hamiltonian is given by 

(9) 

Here cr* {i — x, y, z) are the Pauli matrices and B is the 
spin splitting. The reservoirs {v = L,R) include two 
infinite sets of harmonic oscillators (creation operators 
bi g). Spin-bath interaction strength is denoted by A^^g, 
possibly different at the two ends. 

Let H^j, denotes the local Hamiltonian of the ensemble 
of harmonic oscillators at the u boundary, Hspin = ■j"'^ 
be the Hamiltonian of the spin and Vi^ = \,^,g(T^{bl g + 
bu,q) be the interaction. Using Eq. (8), the energy flux 



from the L contact to the spin unit is given by 

3l = \[T'<^''^^q^LA^)\,^-bL,q) 
q 

+ BayY,^LAbi,q + K<i)l (10) 

or equivalently, 

jL = \[BayXL + a^PL], (11) 

where Xl = Y.q>^L,q{b\^g + bL,q) and Pl = 
i J2q ^L,qOJq{b\ g — bL,q)- An analogous expression exists 
at the R side. It can be shown that the flux operator (10) 
reduces to the stationary heat flux expression utilized in 
Refs. [21, 44, 45] when system-bath couplings are weak 
and the Markovian limit is assumed, 

{jL) = -B[k^^^u-k^^uPd]- (12) 

Here (j) denotes the trace over system and bath degrees 
of freedom, p„ (p^) is the steady state population of the 
up (down) spin level and T^, is the temperature at the v 
contact. The rate constants satisfy the detailed balance 
relation, = k^^j^e~^/'^" , where 

/oo 
e'^-(X.(r)X.(0))rfr 
-00 

= 27r^A2_,[n^(a;,) + l],5(i?-a;,). (13) 
g 

ng{ujg) ~ ^e'^i/T,. _ 2^] ^ jg ^jjg Bose-Einstcin distribu- 
tion function with the Boltzmann constant k-B = 1- 
Equation (12) describes energy current at the L contact 
as the balance between an energy extraction from the L 
reservoir into the spin, and an energy loss from the spin 
to the bath. Appendix B presents in details the deriva- 
tion of this perturbative result from the general operator 
expression (11). 

Similarly, one may analyze the transport properties 
of the diagonally coupled spin-boson model with Vi, = 

Y.q i^i^-A^l.^q + ^t^-.q) ^^'^ Hspin = f o"^ + yc^"^, leading 
to complicated behavior due to the non-separability of 
the two reservoirs [21, 46]. 

Spin-Fermion model. — The spin-fermion model, 
where a spin impurity is coupled to two Fermi seas of 
different temperatures and/or chemical potentials, is an- 
other example of a hybrid structure, useful e.g. for study- 
ing electronic and radiative heat transfer between metals 
[22], 

HsF = -7^a'- ^^ekc\^^Cy^k + X! "''.fc.94,feCi.,g-(14) 

The first term here accounts for spin splitting. The sec- 
ond term includes the two independent reservoirs (leads) 
of spinless electrons, creation operator c|,^, iv = L,R). 
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We assume that the leads are kept (each) in thermal equi- 
librium at temperature Tj^ and chemical potential ii^. 
The last term in (14) describes spin-bath interactions, 
where we disregard charge tunneling between the metals 
and allow only for transfer of energy excitations. Utiliz- 
ing Eq. (8), the heat current at the L contact is given 
by 



k,q 



(15) 



If the metals have strictly linear dispersion relation this 
result can be exactly mapped into a bosonized descrip- 
tion [48] to yield the current (10). Deviations are ex- 
pected when the metals have energy dependent density 
of states [22]. Following the derivation sketched in Ap- 
pendix B, taking into account the fermionic nature of the 
operators, one can show that in second order system-bath 
coupling, going into the Markovian limit, the stationary 
heat current is given by Eq. (12) with the rates 

/oo 
e-^-(F.(r)F.(0))dT 
-OO 

= 27r^ |a^,fe,g|^n^(efc)[l - n^(e,)]5(efe - e, - B) 

k,q 

kUu = K^ae-'"^^, (16) 

where = J2k,q ^f.fe.^^t.fe^i'.?' force the bath ex- 

erts on the system, and np{e) = [e(^~^'^)/^" -|- 1]"-*^ is the 
Fermi-Dirac distribution function of the ly bath. 

The perturbative rate expression (12) also holds for 
mixed boson-fermion systems, e.g. when energy is di- 
rected from a phonon bath into an electronic excitation 
through a local impurity. One simply employs then the 
expressions (13) and (16) for the phononic and electronic 
bath-induced transitions. 



show that ji^^ = 0, thus the flux operator is given by 



Js = 



dV{s,s + l) 



dV{s,s + l) 



dx 



s+l 



(17) 



where {}_,_ denotes the anticommutation relation. This is 

just the quantized form of the classical flux defined in ref- 
erence [25]. For the the quadratic interaction {xs—Xs+i)"^ 



we should exclude the local terms xi and 



shift 



them into /i° and respectively, as discussed in Sec- 
tion II. For a bilinear coupling model we thus consider the 
interaction V{s, s + 1) = XxgXs+i with spring constant 
A. The flux operator then reads 



js = aX{psXs+i - XsPs+i)/'2 
= -iaX{blbs+i - &^+i6a)/2, 



(18) 



where the second line is the bosonic expression with the 
creation (annihilation) operator 6| {bg). The commu- 
tation relation between the total Hamiltonian and the 
current operator is given by 



Xs+l- 



dUjx, 
dx. 



0(A2 



Therefore, if 



dUjx) _ 
dx 



(19) 

= 0, or oc x^, [js,H] = within 
first order coupling. This implies that free particle mo- 
tion and harmonic potentials pertain a constant current, 
or in other words, the heat current is conserved in these 
systems. 

One can also calculate the higher order term in (19), 
O(A^) = iaX'^{xl + XsXs+2 - Xs-iXg+i - xl_^_i)/2. If the 
total flux is defined as J = J2'^=i isi the commutation re- 
lation between the complete flux and the total Hamilto- 
nian is given by [J, H] — iX^ | {xf — x^) for the harmonic 
U oc .T^ potential. Therefore, the Heisenberg equation of 
motion reads 



aX^ 

dJ/dt = — -— (xf — x%). 



(20) 



IV. CURRENT CONSERVATION CONDITIONS 
FOR BOSONIC AND FERMIONIC SYSTEMS 



With the help of the heat flux operator we can obtain 
general properties of specific quantum systems [39]. This 
is in contrast to standard calculations where one needs 

to make use of specific quantum states [33-35] . We prove 
next that linear harmonic systems and some special spin 
chains (XY, Ising) have zero thermal resistance using the 
operator form of the energy flux. 

Bosons. — We consider the quantized system used in 
[25], /ig = p'^/2 + U{xs) and inter-site potential V{s, s + 
1) = V{xs, Xg^i), where Xs and ps are the coordinate 
and momentum of the particle at the s site. It is easy to 



This result shows that the total flux depends only on the 
contacts properties: coupling strength and temperature 
(going into thermal averages). Furthermore, in closed 
loop systems, the complete current J is a constant op- 
erator. This conclusion is well established, however, we 
give here a simple proof of the operator form, without 
the need to go into the system's quantum states. It 
can be shown that the current is also conserved for dis- 
ordered ID harmonic systems. For example, assuming 
different force constants between sites As^s+i, one gets 
dJ/dt = f (Af 22^1 - ^'n-i.n^'n)- 

Fermions: Nearest Neighbor Spin systems. We con- 
sider next a periodic spin chain of length A''. The sys- 
tem can be mapped into a system of fermions using the 
Wigner- Jordan transformation, see e.g. [49] . Let the on- 
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site potential /i° and the inter-site potential V be 

K - ifs 

V{s, s + 1) = A(A<7>,\i + Ba\al^^ + Calal^^\ 

(21) 

where A, B, C arc the interaction coefficients. It is easy 
to show that the first-order flux operator is given by 



,•(2) 

Js-ys+1 



aeX 



A + B 



(22) 



Using the Wigner- Jordan transformation, the current can 
be also rewritten as 

= iaeX{A + B){cl^,c, - clc,+,), (23) 

expressed in terms of spinless fermionic creation and an- 
nihilation operators c| and Cg respectively. The second 

order contribution js'^^ oc is nonzero in general, but is 
too cumbersome to be included here. 

(2) 

The current operator jg_^s^i is essentially the stan- 
dard spin current operator multiplied by the bias e. This 
term reflects energy flow due to spin current, while ji^' 
accounts for thermal energy flow [39]. At weak inter- 

(2) 

site coupling, A -C e, dominates the energy cur- 

rent, while for zero magnetic fields only ji*^ survives. 
Throughout the paper we always assume nonzero mag- 
netic splitting e, imlcss otherwise stated. 

We continue and analyze current conservation in the 
model (21), 

[js,H] = /i° + + n^Kh!^ + h"s+i] 

+ [j?ls+i,V{s -l,s) + V{s, s + l) + V{s + l,s + 2)]. 

(24) 

To the first order in A the commutator is therefore given 

by 

[-,•(2) ftO -L ftO 1 

lJs^s+1 ' "-s + "-s-l-lj 

= ae'X^^[al^,a^, - ay^^.+^al + al^,] = 0. 

(25) 

Thus, for the periodic spin chains considered here, only 
high order terms in A may lead to current decay. We dis- 
cuss next some special cases: (i) A = —B, corresponding 



,(2) 



0, im- 



to the antifcrromagnctic phase. Here jl^^^^ 
plying that there is no current in the antiferromagnetic 
phase in the first order approximation, (ii) The Heisen- 
berg model, A = B = C . In this case the flux operator 
(6) agrees with the definition of Ref. [29], see also Ap- 
pendix A, since [{a^+ag_^_i,V{s,s + l)] = 0. This system 
was extensively investigated in Refs. [29, 38]. (iii) The 
XY model, A = B and C = 0. We calculate here the 
high order contribution to the current and find 

.aAM2 



.•(4) 



-(4c 



s+2 



+cJ_iCs+i - 4+iCs_i), 



(26) 



Combining Eq. (23) with Eq. (26) we get that [js,H] = 
+ O(A^) in the XY model. The current operator is 
therefore a constant in the first order approximation, 
while the total current exactly becomes 



dJ/dt = 8aeA^A^(ni — njv), 



(27) 



in analogy with Eq. (20) for the bosonic Hamiltonian. 
Here rig = cjcg is the number operator. We conclude 
that the total current across the systems depends only 
on the properties of the chain's ends. Thus, in closed 
loop systems the total current J is a constant operator. 
As a final (iv) case we consider the transverse Ising 



model, B 



C = 0. 



cfcTs+i), js^'' = 0. The commutator [js, H] is zero in first 
order of A while the second order contribution, resulting 

from the commutator [^g^-s-i-i > ~ 1> s) + V{s, s -|- 1) -|- 
V{s + l,s + 2)], leads to 



Here 



dJ/dt = 4aeA^A^(ni - njv). 



(28) 



We can summarize our observations as follows: If a 

Hamiltonian is written by a linear combination of bilinear 

operators, a bosonic set |fos&t, bib], b^bt^ and a fermionic 

set |cjci, c^cj, c,5Ct|, it can always be expressed in terms 

of quasiparticle operators jq, where H = X^^^ ^qlllq 
(see e.g. [50]). Since there is no interaction between the 
quasiparticlcs, the systems behaves like a collection of 
free particles. The harmonic oscillator chain with linear 
couplings is an example of bosonic Hamiltonian. The 
XY models are examples for independent fermions. Both 
systems yield ballistic motion with no thermal resistance. 
In contrast, the Heiscnberg model with nonzero magnetic 
field does not belong to such systems because it contains 
an on-site interaction clcgclcf when C is not zero [20]. 



V. NECESSARY CONDITION FOR CURRENT 
CONSERVATION FOR A GENERAL 
ONE-DIMENSIONAL SYSTEM 

We Consider a chain of length TV with M levels at 
each site. The commutation relation between the total 
Hamiltonian and the flux operator can be written as 



[js,H] = FiX) + OiX'), 



(29) 



where F{X) is the first order term and O(A^) contains 
higher order terms in A. The necessary condition for 

current conservation is F(X) = 0. Wc emphasize that 
this is only a necessary condition. If F(X) ^ the system 
potentially shows a diffusive dynamics. 

The most general Hamiltonian for this system can be 
generated by generator set g = { h , E""}, h = 
(n^, n^, n^) is the vector operator with n' = \i) 
E" denotes M"^ — M operators |z) {j\ where i j. The 
vectors ~a^s are M-dimensional, and are usually referred 
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to as roots [51]. The commutation relation between h 
and E°' is [h ,E"] = 'aE " , where the vector 7? can 
be considered as an eigenvalue of the vector operator 
h. For instance, in the two level system (M = 2), 
t = (|1)(1|,|2)(2|), E-^ = (|1)(2|,|2)(1|) and there 
are two roots ~ai — (1,-1) and 0:2 = (—1)1)) corre- 
sponding to E'^^ = |1) (2| and = |2) Appendix 
C presents the M = 3 case. 

Using this notation, the most general Hamiltonian up 
to a two-body interaction can be written as 



(30) 



where the vector e 



energy level and V-^ 



(ei, £2, em), is the i state 
are inter-site coupling parame- 
ters. The units are assumed to have identical spectra 
and wc use constant nearest-neighbor interactions along 
the chain. The last term in (30) includes many body in- 
teractions Vd{s, s + 1) = ^ Ui^jn\rv'g_^-^. It is easy to show 
that the commutator of the current with H yields 



3+V 



(31) 



The necessary condition for current conservation, F{X) 
0, therefore implies 



±~e ■ (3 



a, /3 ■ 



(32) 

This condition 



for nonzero coupling parameters V- 

(with the plus sign) is naturally fulfilled for harmonic 
systems, since (ej_i — ej) — {ek-i ~ ffc) for any j, k. For 
fermionic models M=2 and the • /3 = —~e ■ ~a condi- 
tion is trivially conformed. Both systems indeed lead to 
current conservation, see Section IV. 

For a system with an arbitrary spectra this condition 
translates into "a = ± /? , implying that the interaction 
contains only the following terms: E^ E^_^^, E!^ EJ^^ 
and Vd{s, s + 1) [52]. The corresponding current operator 
is 



-1-1 • 



(33) 



This expression reduces into the fermionic limit (Section 
IV) when M = 2. The M = 3 case is exemplified in 
Appendix C. 

The necessary condition (32) is an imperative step to- 
wards identifying normal transport (Fourier) systems, as 
it helps us pinpoint current conserved systems directly, 
without detailed numerical calculations. If the system 
satisfies • c? 7^ zb"? • 13 , one can directly deduce that 
the thermal current is not conserved. Note that in the 
Heisenberg model F{X) = 0, and only the next term in 
Eq. (29) is finite, accounting for dissipation of energy 
[53]. 



VI. FORMAL FOURIER'S LAW 

Recently, there are several ideas of how to approach 
Fourier's law from fimdamcntal principles [26- 30]. Here 
we will show that the appropriately defined flux opera- 
tor naturally leads to the discrete form of the law. The 
derivation yields the conductivity coefficient for a general 
ID system in terms of the Hamiltonian parameters. We 
begin with a generic nearest-neighbor Hamiltonian 



(34) 



including local interactions and inter-site couplings. In 

our definition (6), the average flux, j = Tr{pj}, at weak 
interactions reads 



7(2) 



la 



Tr{p[{h°+i-h%V{s,s + l)]} 



-^Tr{Ah°T{t)}, 



(35) 



using the cyclic property of the trace. Here p is the total 
density matrix, A/i° = — /ig is the difference between 
local energies at neighboring sites and T{t) = i[V{s,s + 
l),p(t)] is hermitian. We can also write this expression 
explicitly in terms of local s functions. 



s+1 



,)) 



(36) 



where = Tr{h'^T(t)}. If wc define a local temperature 
Tg at each site, we can then relate the current between 
sites with the temperature difference ATg = Tg+i — Tg, 



7(2) 



= - a 



At 



Cs 



An 



(37) 



where Ag^ = 5^+1 - g^, and 



A7^ 

AT, 



is the specific 

heat. This is the discrete Foiuier's law [29, 38]. We 
can identify the microscopic-local thermal conductivity 
as Kg ^— Qi -;^Cs, as long as g^ can be uniquely defined 



(see discussion below), and the ratio is finite. 

As an example we consider a three-spin system. For 
the XY model, if the initial state is lO)! ll^olC),, it is 

2V2A 



easy to show that 



/I 



A/12 
A/i, 



2 l"/3' 

. For weak 



3cos2y2A«+l 

2\^t holds. The heat conduc- 



couphng. At < 1, —ik 

^ ' a4 

tivity is then given by k = 2X^tCs, in agreement with 
our recent calculation [30]. It also shows that although 
the total current of the XY model is conserved, the par- 
tial current between two sites may have the form of the 
Fourier's law before thermal eqiulibrium sets [24]. 

We explain next how to define g^ uniquely. Although 
we could formally write Eq. (37) , g^ may not be uniquely 
defined because T{t) depends on the index s: T{t) could 
be either defined as ^[^^(s, s+l),p{t)] or z[F(s, s—l),p{t)]. 
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Therefore, the condition for to be exclusively defined 
is 



APPENDIX A: An Alternative, symmetry limited, 
definition for tlie lieat current 



TV,{[V^(s,s + I), Pirn = Tr,{[V{s, s - (38) 

The trace Tr,, runs over all sites except site s. It is easy 
to show that Ps-i,s+ip{t)Ps-i,s+i = p{t) is a sufficient 
condition for satisfying Eq. (38), where Pg-i^g+i is the 
exchange operator between sites s — 1 and s + 1. If the 
total Hamiltonian is invariant under Pg-i^s^i, as it is in 
many physical cases, the last condition translates into a 
condition on the system preparation. 



-l,s+l 



s-l,s+l 



p(0). 



(39) 



This is a sufficient (but not necessary) condition for at- 
taining a unique expression for g^. Once g^ is carefully 
defined, we can proceed and calculate the thermal con- 
ductivity using Eq. (37). In the example above the initial 
state was set to |0)j^ 1 1)2 10)3, which is indeed invariant 
under the exchange Pi, 3. Note that since the validity 
of the Fourier's law is independent of initial conditions, 
the requirement to fulfill Eq. (39) is solely meant for 
distinctively identifying the conductivity. 



VII. CONCLUSION 

In this paper we present and re-examine the heat flux 
operator that exactly satisfies the continuity equation 
for a general Hamiltonian in one dimension. Based on 
the definition, we deduce the necessary conditions on the 
inter-site interaction that result in current conservation. 
This analysis sets the first step towards the exploration of 
the validity of Fourier's law of heat conduction in Hamil- 
tonian systems: systems that conserve energy have di- 
verging conductivity. As an example, using a simple op- 
erator algebra, we prove that independent bosons and 
fermions conduct heat ballistically. We further apply 
the definition to various impurity models, relevant for 
understanding heat flow in nanojunctions, and obtain a 
non-perturbative non-stationary expression for the heat 
current. The microscopic heat conductivity coefficient 
naturally emerges in the present definition. 

While previous works have typically relied on specific 
quantum states, calculating only expectation values, see 
for example [34-38] , the results presented here essentially 
depend only on operator calciflations. Possible exten- 
sions include generalization of the heat current definition 
to time dependent situations, and exploration of the nec- 
essary condition for the applicability of the Fourier's law 
of heat conduction in ID chains [30] . 
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We follow here a symmetry-limited definition of the 
heat flux operator often adopted in the literature [29, 36, 
37] . The generic ID Hamiltonian includes local potentials 
and inter-site interactions 



(Al) 



The heat flux operator is defined by considering the time 
evolution of the local, non-interacting energy operator, 

= -iK V{s, » - 1)] - ilh», V{,, s + 1)]. (A2) 

We next assume that a continuity equation for holds, 
based on the approximation that the local energy is con- 
served [36] 



dh^s ^ Us-i ~js) 
dt a 



(A3) 



By comparing Eq. (A2) with (A3) one can identify the 
current between sites as 



ia[hlV{s,s + l)]; 



Js 



-ia[hlV{s,s-l)]. 



(A4) 



However, the second equality above produces js = 
—ia[h^g^i,V {s + l,s)] when shifted to site s. This can 
be consistent with the first equality of Eq. (A4) only if 
the condition 



1)1-0 



(A5) 



is satisfied. The definition (A4) is thus restricted to a 
limited class of Hamiltonians that satisfy (A5). We em- 
phasize again that the heat current was defined here by 
studying local, non-interacting energy changes, while Eq. 
(5) defines the heat current by studying the total en- 
ergy at a site, incorporating intcr-sitc interactions, see 
Eq. (4). The Heisenberg spin-^ model, /i° = |crf, 
V{s,s + 1) = A(a>f+i -I- (jyal^i + o-^+i), is an ex- 
ample of a system obeying (A5). 



APPENDIX B: Spin-Boson model: Derivation of 
the weak coupling expression for the heat current 

We derive here a weak-coupling expression for the 
steady-state heat current in the spin-boson model us- 
ing the non-perturbative definition (10). The two-bath 
(y = L, R) spin-boson Hamiltonian is given by 



HsB = Hgpin -\- ^ -I- ^ K, 



(Bl) 
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where 



H. 



B 



spin 



a'; Hl^ = J2^gbl<ibu,g; = a'^X^. (B2) 



Here B is the spin spHtting, bl ^ is a creation operator 
satisfying the bosonic statistics, and Vi, includes system- 
bath interactions at each contact, X^, = K.q{bl_q + 
bu.q)- There is no direct coupling between the two har- 
monic baths (temperature T^,), as they are coupled only 
through the central spin. 

The general expression for the current operator is given 
by Eq. (6), jl = ^Wl - iJ^pm,^!,], disregarding for 
convenience the factor a. Note that j^^^ = 0, see Eq. 
(7), since [Vl,^/?] = 0. In the present model the current 
operator from the L interface to the spin is given by 



where Pl = i^q'^L,q^q{b\ ^ ~ bL,q) denotes the sum 
of the momenta of the harmonic oscillators at the 
left boundary. This expression is valid in the non- 
perturbative regime and for non-stationary situations. In 
steady-state the expectation value of the interaction is 
zero, e.g. at the L contact. 



(B3) 



OVl 
dt 



(B4) 



Since = —Ba^ and Xl = Pl, find that {(J^Pl) = 
{BcjyXi). The stationary heat current is therefore given 

by 



{jl) = Tr{pji} = BTvipayXL}, 



(B5) 



where p is the total density matrix. Using the energy 

representation, = \u){u\ — \d){d\, = \d){u\ + \v,){d\, 
= —i\u){d\ +i\d){u\, we can write the heat current as 



Oi) = iBTrB{{Pu,d - Pd,u)XL}, 



(B6) 



where Tre denotes the trace over the thermal baths (L 
and R) states only. This expression can be evaluated by 
solving the Liouville equation, written here explicitly for 
the nondiagonal matrix element 

Pd,u{t) = iBpd,u{t) - iX{t)pu,u{t) + ipdAt)X{t), (B7) 

withX = Xl+Xr. Formal integration of this differential 
equation yields 

PdAt) = f e^^(*-")[-iX(r)p„,„(T) + ipd,d{T)X{T)]dT. 
Jo 

(B8) 

We evaluate next the term TtbIpci^uXl} under the fol- 
lowing approximations: (i) weak system-bath coupling, 

neglecting higher order correlation functions, (ii) Marko- 
vian limit, assuming the spin's relaxation timescale is 



longer than that of the bath fluctuations, and (iii) ini- 
tial factorized condition, where p is well approximated 
by the product p{t = 0) = Pspin{t = 0)plPr. Here 
p^ ~ e~^'^/-^'^/Tr{e^^''/-^"} are the density operators of 
the thermal baths. These assumptions are compatible 
with the Redfield approximation [54]. Using (B8) we get 



TrB{pd,«XL} 



+ Wd 



iPu{t) / e*^"(Xi(T)Xi(0))dT 
Jo 

oo 

Bt I 



{t) r 

Jo 



e'^^(Xi(0)Xi(r))dr, 



(B9) 

where p„ = Ttb{Pu,u} denotes the population of the spin- 
up state and pd is the spin-down population. Note that 
terms of the form (Xi(t)X/{(T)) disappear, since the two 
reservoirs are not correlated. Following the same proce- 
dure for the second term in Eq. (B6) we obtain 



TrB{p„,dXi,} = ipu{t) / 

J —OC 

iit) f 

J —OC 



lpd{ 



e'^"(Xi(T)Xi(0))dT 



^"^XL{0)XL{T))dT. 



(BIO) 



Combining equations (B9) and (BIO) provides us with 
the stationary thermal current under weak-coupling and 
Markovian approximations, 

(ii) = -B\puk!^^d-Pdkd^uh 
with the relaxation rates 



-{XL{T)XL{0))dT 



(Bll) 



OO 
OO 



k: 



e-^^-(X.(T)X.(0))dT. (B12) 

-OO 



Equation (Bll) describes energy current through the 
jmiction, calculated e.g. at the L contact, as the bal- 
ance between an energy gain from the reservoir to the 
spin, and an energy loss from the spin to the L bath. 

The diagonal elements of the density matrix, pd and 
Pu, can be further calculated under the same set of ap- 
proximations, to yield the quantum Master equation, 

Pu = -Pu{t) k'^^^d +Pd{t) XI ^d^u 

Puit)+Pd{t) = l. (B13) 
In steady state {p = 0) the spin occupations are 



Pu = 



Pu +Pd = 1- 



k? 



h.R 



UR ' 
'^u^d 



(B14) 

Plugging Eq. (B14) into (B12) leads to an explicit ex- 
pression for the current 

\Jl) - tr — —. r-TR — —. — • (^15) 



k^-^d ■ 



f^u^d 



k^_ 



^d^u 
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An analogous expression holds at the R contact. This 
is the well established quantum Master-equation limit, 
used in various apphcations [21, 22, 44, 45, 55-57]. 

We can also extend the calculations to non-stationary 
situations. In this case one needs to evaluate the extra 
term Trja^Pi,} = TrB{(p„,d + Pd,u)PL} in Eq. (B3), 
resulting in momentum-position correlation functions of 
the form, {PL{t)XL{T)) in second order system-bath cou- 
plings. 

Note that [H^ + i?spi„, Vl] ^ for the spin-boson 
Hamiltonian. Therefore, we cannot use in general the 
definition of Appendix A, = *[-ffspMi, Vl]. This lim- 
ited expression is still applicable in a steady-state situa- 
tion since {dVL/dt) = 0, translating into {[Hsb,Vl]) = 
{[Hi + Hspin, Vl]) = 0, see Eq. (A5). 

APPENDIX C: Current conservation in an M=3 
states model 

We clarify the notation and the results of Section V 
using an M=3 level system. According to our notation, 
the diagonal operators are 

7r = (|l)(l|,|2)(2|,|3)(3|). (CI) 

The six nondiagonal operators with their respective roots 
are 





= |1) (2|, 


"ai 


= (1,-1,0) 




= |2) 


^2 


= (-1,1,0) 




= |1) (3|, 


"as 


= (1,0,-1) 


E"^ 


= |3)(1|, 




= (-1,0,1) 


E"^ 


= |2) (3|, 


"as 


= (0,1,-1) 




= |3) (2|, 




= (0,-1,1) 



(C2) 



where the energies at each site are ~€ = (ei,e2,e3). 
If the system conserves current [i.e. it fulfills (32)], 
the site-site interaction can include only the following 
terms: Ef^Ef_^_\, (n = 1..6), and the pairs Ef^Ef^j^, 

Ef'> Ef±^\ , Ef"- Ef^-^ . The current operator in this model 
is given by Eq. (33), 

= ia|^^„^,(ei - e2){ErEf_^\ - Ef^Ef_^\) 
+V^,,1SS^2 - e3)(i^r i^Ji - Ef^Ef_^\)Y 

(C3) 

a generalization of the spin chain result (22). 
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